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ABSTRACT 



Aims. This article shows the first evidence for gravitational lensing phenomena in high energy gamma-rays. This evidence comes 
from the observation of a gravitational lens induced echo in the light curve of the distant blazar PKS 1830-211. 
Methods. Traditional methods for the estimation of time delays in gravitational lensing systems rely on the cross-correlation of the 
light curves of the individual images. In this paper, we use 300 MeV-30 GeV photons detected by the Fermi-LAT instrument. The 
Fermi-LAT instrument cannot separate the images of known lenses. The observed light curve is thus the superposition of individual 
image light curves. The Fermi-LAT instrument has the advantage of providing long, evenly spaced, time series. In addition, the photon 
noise level is very low. This allows to use directly Fourier transform methods. 

Results. A time delay between the two compact images of PKS 1830-21 1 has been searched for both by the autocorrelation method 
and the "double power spectrum" method. The double power spectrum shows a 3 a evidence for a time delay of 27.5+1.3 days, 
consistent with the result from Lovell et al. (1998). The relative uncertainty on the time delay estimation is reduced from 20% to 5%. 

Key words. Gravitational lensing: strong - [Galaxies] quasars: individual: PKS 1830-21 1 - Methods: data analysis 



1. Introduction 

The precise estimation of the time delay between components 
of lensed Active Galactic Nuclei (AGN) is crucial for model- 
ing the lensing objects. In turn, more accurate lens models give 
better constraints on the Hubble constant. More than 200 strong 
lens systems have been found. Most of them were discovered 
in recent years by dedicated surveys such as the Cos mic Lens 
All-Sky Survey dMvers et alJl2003b iBrowne et alJl2003l) and the 
Sloan Lens AC S Survey dBolton et al.ll2004l) . The launch of the 
Fermi sattelite (lAtwood et al.ll2009h in 2008 gives the opportu- 
nity to observe gravitational lensing phenomena with high en- 
ergy gamma rays. The observation strategy of Fermi-LAT, which 
surveys the whole sky in 190 minutes, allows a regular sampling 
of quasar light curves with a period of a few hours. Since the 
launch of the Fermi satellite, the LAT instrument has been col- 
lecting high energy photons for more than 800 days. This pa- 
per deals with the observation and estimation of time delays of 
strong lensing systems and not with the detectio n and use of mi- 
crolensing phenomena as in Torres et al] d2003h . 

The multiple images of a gravitational lensed AGN cannot be 
directly observed with the high energy gamma-ray instruments 
such as the Fermi-LAT, Swift or ground based Cerenkov tele- 
scopes, due to their limited angular resolutions. The angular res- 
olution of these instruments is at best of a few arcminutes (in 
the case of HESS), when the typical separation of the images for 
quasar lensed by galaxies is a few arcseconds. Since the multiple 
images cannot be observed, methods based on the detection of 
an echo in the light curve are studied in this paper. The observed 
signal is the superposition of the intrinsic light curve of the AGN 
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and at least one delayed copy of this light curve, with a different 
amplitude. 

The method has been tested on simulated light curves and on 
Fermi observations of the very bright radio quasar PKS 1830- 
211, for which the time delay has been accurately measured by 
lLovell et"aT1 (Q998) using radio observations. The paper is orga- 
nized as follows: we first give a very brief summary of proper- 
ties of PKS 1 830-2 1 1 and of Fermi LAT data towards this AGN. 
Then we introduce the methods for time delay estimation. The 
last section is devoted to the measurement of the time delay be- 
tween the two compact components of PKS 1830-21 1. 



2. The PKS 1830-211 gravitational lens 

The AGN PKS 1830-211 is a variable, bright radio source 
and an X-ray blaza r. Its redshift was measured to be z=2.507 
dLidman et al.f l999). The blazar was detected in the y-ray wave- 
lengths with EGRET. The ass ociation of the EGRET source with 
the radio source was done by iMattox et al.l dl997l) . The classifi- 
cation of PKS 1830-211 as a g ravitational lensed quasi-stellar 
object was first proposed by iPramesh "Ra o & Subrahm anvanl 
dT988l) . 

Two potential lenses have been identified along the line of 
sight to PKS 1830-211 by the observation of molecular ab- 
sorption at centimeter and millimeter wavele ngths. One possibl e 
lensing object is a galaxy at redshift z=0.19 (Lovell et al. 1996), 
and t he other a gala xy at redshi f t z=0 .89 (Wiklind & Combes 
1996). According to IWinn et all d2002l) . the actual lens is the 
galaxy at redshift z=0.89. 

PKS 1830-211 is observed in radio as an elliptical ring- 
like structu re connecting 2 bri ght sources distant of roughly one 
arcsecond dJauncev et al.ll 199 lb . The compact components were 
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Fig. 1. Fermi LAT light curve of PKS 1830-21 1, with a 2 days binning. The energy range is 300 MeV to 300 GeV. 



separately observed by the Australia Telescope Compact Array 
at 8.6 GHz for 18 mont hs. These observ ations and the subse- 
quent analysis made bv lLovell et al.l (1199 8) gave a magnification 
ratio between the 2 images equal to 1 .52 + 0.05 and a time delay 



The power spectrum P v of the source is obtained by comput- 
ing the square modulus of g(v). 



of 26^ days. 



2. 1. Fermi LAT data on PKS 1830-21 1 

PKS 1830-21 1 has been detected by the Fermi-LAT instrument 
with a detection significa nce above 41 Ferm i Test Statistic (TS), 
equivalent to a 6 <x effect ( Ab do et al.l 2010). The long-term light 
curve is presented on figure Q] with a 2 days binning. The data 
analysis from this paper uses a 2 days binning, which provides a 
sufficient photon statistic per bin with a time span per bin much 
shorter than 28 days. The data analysis was cross-checked by 
binning the light-curve into 1 day and 23 hours bins, with sim- 
ilar results. The data were taken between August 4 2008 and 
October 13 2010, and processed by the publicly available Fermi 
Science Tools version 9. The v9rl5p2 software version and the 
P6_V3 -DIFFUSE instrument response functions have been used. 
The light curve has been produced by an aperture photometry 
selecting photons from a region with radius 0.5 deg around the 
nominal position of PKS 1830-211 and energies between 300 
MeV and 300 GeV. 



3. Data processing and method 



(v)| 2 



|/(v)| 2 (l + b 2 + 2bcos(2nva)) 



(1) 



The measurable P v is the product of the "true" power spec- 
trum of the source times a periodic component with a period (in 
the frequency domain) equal to the inverse of the relative time 
delay a. 

The usual way of measuring the time delay a is to calcu- 
late the au tocorrelation function of f( t). This method was inves- 
tigated by Gei ger & Schneiderl(ll996l) . The autocorrelation func- 
tion can be written as the sum of three terms. One of this term 
describes the intrinsic autocorrelation of the AGN, decreasing 
with a time constant A. If A is larger than a, the autocorrelation 
method fails, because the time delay peak merges with the in- 
trinsic component of the AGN. Another potential problem with 
the autocorrelation method is the sensitivity to spurious period- 
icities such as the one coming from the motion and rotation of 
the Fermi satellite. 

The periodic modulation of P v suggests the use of another 
method, based on the computation of the power spectrum of P v , 
noted D a . This method is similar in spirit to the cepstrum analy- 
sis dBogert et alJI 1963b used in seismology and speech process- 
ing. If |/(v)| 2 was a constant function of v, D a would have a 
peak at the time delay a. In the general case, D„ is obtained by 
the convolution of a Dirac function, coming from the cosinus 
modulation, by the Fourier transform of the function: 



3.1. Idea 



h(v) = f(v)Window(0, W) 



(2) 



If a distant source (in our case an AGN) is gravitationally lensed, 
the light reaches the observer through at least 2 different paths. 
We assume here that there are only 2 light paths. In reality, the 
light curves of the 2 images are not totally identical since (in 
addition to differences due to photon noise) the source can be 
microlensed in one of the two paths. In this paper, we neglect 
this effect. Neglecting for the moment the background light, the 
observed flux can be decomposed into 2 components. The first 
component is the intrinsic AGN light curve, given by /(f), with 
Fourier transform f(v). The second component has a similar 
time evolution than the first one, but is shifted in time with a 
delay a. In addition, the brightness of the second component dif- 
fers by a factor of b from that of the first component, so that it 
can be written as bf(t + a) and its transform to the Fourier space 
gives bf(v)e- 2niva . 

The sum of two component gives g( f) - f( f)+bf(t+a) which 
transforms into g(v) = /(v)(l + be~ 2mva ) in Fourier space. 



where Window(0,W) is the window in frequency of P v and 
W is maximum available frequency. The Fourier transform h(a) 
of h(v) defines the width of the time delay peak in the double 
power spectrum D„. 

For instance if /(v) = e ~ lv and AW » 1 , then the time delay 
peak in D a has a Lorentzian shape with a FWHM of A/n. For a 
typical value of A — 10 days, one has a FWHM of 3 days. 

In next section we describe the calculation of P v and D a and 
illustrate the procedure with Monte Carlo simulations. 

3.2. Power Spectrum 

The Fourier transform is a powerful technique to analyze astro- 
nomical data, but it requires a proper preparation of the observa- 
tions. To avoid problems arising from the finite length of mea- 
surements, sampling a nd aliasing, we use the p rocedure for data 
reduction described bv lBrault & Whitd dl97lh . 
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Fig. 2. Simulated power spectrum of a lensed AGN with a time 
delay between images of 28 days. The power P v is in arbitrary 
units. 



We start with dividing the whole light curve into several seg- 
ments of equal length. The data have to be evenly spaced and 
the number of points per segment needs to be equal to a power 
of 2. The choice of the segment length is compromise between 
the spectral resolution and the size of error bars on points in the 
spectrum. The resolution of lines increases with the number of 
points in the segment, but the error on each point in the spectrum 
decreases as the number of segments. 

As suggested by iBrault & White! (Il97lh . the segments are 
overlapped to obtain a larger number of segments with a suffi- 
cient number of points. Then the data in each segment are trans- 
formed with the following procedure: 

1 . Data gaps are removed by linear interpolation 

2. The mean is subtracted from the series to avoid having a 
large value in the first bin of the transform 

3. The data are oversampled to remove aliasing. 

4. Zeros are added to the end of the series to reduce power at 
high frequencies 

Finally, the power spectrum is averaged over all segments. 
The power spectrum obtained with an artificial light curve is 
shown on figure [2] 

The artificial light curve was produced by summing three 
simulated components. The light curve of PKS 1830-21 1 shown 
on figure[T]does not exhibit any easily recognizable features, but 
has a rather random-like aspect. The first component was thus 
simulated as a white noise, with a Poisson distribution. It would 
be more realistic to use red noise instead of white noise but the 
latter is sufficient for most of our purposes, such as computing 
D a . The second component is obtained from the first by shifting 
the light curve with a 28 days time lag. The effect of differential 
magnification of the images has also been included. The back- 
ground photon noise was taken into account by adding a third 
component with a Poisson distribution. 

The mean number of counts per 2 day bin for PKS 1830- 
21 1 is 5.42 counts. This value was used for the simulation of the 
artificial light curve. The first and second component contribute 
80% of the simulated count rate and the rest is contributed by 
the noise component. 

To get a simulated flux similar to the observed flux, the sim- 
ulated count rate is divided by the average exposure from obser- 
vations ( 2.8557 1 7 s cm 2 ). 



3.3. Time Delay Determination 

The methods of time delay determination use the power spec- 
trum P v obtained as described in the previous section. The sim- 
ulated P v presented on figure [2] shows a very clear periodic pat- 
tern. From equation Q] we know that the period of the observed 
wobbles corresponds to inverse of the time delay between the 
images. 

Our preferred approach was to calculate the double power 
spectrum D a . As in section 13.21 the power spectrum P v has to 
be prepared before undergoing a Fourier transform to the "time 
delay" domain. 

The low frequency part (v < 1 /55day~') of P v is cut off. This 
cut arises because of the large power observed low frequencies in 
the power spectrum of PKS 1830-21 1. The high frequency part 
of the spectrum P v is also removed because the power at high 
frequency is small (it goes to at the Nyquist frequency). The 
calculation of D a proceeds like in section [3~2l except that the P v 
data are bend to zero by multiplication with a cosine bell. This 
eliminates spurious high frequencies, when zeros are added to 
the P v series. The D a distribution is estimated from 5 segments 
of the light curve. In every bin of the D a distribution, the esti- 
mated double power spectrum is given by the average over the 5 
segments. The errors bars on D a are estimated from the disper- 
sion of bin values divided by 2 (since there are 5 segments). Due 
to the random nature of the sampling process, some of the error 
bars obtained are much smaller than the typical dispersion in the 
D a points. To take this into account, a small systematic error bar 
was added quadratically to all points. The result (with statistical 
error bars only) is presented on figure[3] 

As described in section [3~2l we simulated light curves with 
a time delay of 28 days. A peak is apparent near a time delay of 
28 days on the D a distribution shown on figure[3] The points just 
outside the peak are compatible with a flat distribution. Including 
also the points in the peak gives a distribution which is incom- 
patible with a flat distribution at the 10 sigma level. 

The parameters of the peak were determined by fitting the 
sum of a linear function for the background plus a Gaussian 
function for the signal. In the case shown on figure [3] the time 
delay estimated from D a is 27.94 + 0.61 days. 

As mentioned in section [3~T1 the usual approach for the time 
delay estimation is to compute the autocorrelation of the light 
curve. The auto-covariance is obtained by taking the real part 
of the inverse Fourier transform of P v . The auto-covariance is 
normalized (divided by the value at zero time lag) to get the au- 
tocorrelation. The autocorrelation function of an artificial light 
curve simulated as in section [3~2l is presented on figure[4] 

A peak with a significance of roughly 16 cr is present at 
27.85 + 0. 14 days. However the significance of this peak is over- 
estimated since light curves are simulated with white noise in- 
stead of red noise. The autocorrelation function of a light curve 
driven by red noise is given by e~ fl//1 . In the case of our simu- 
lated light curves, A — 0, so that the peak is a little affected by 
the background of the AGN. 

For the simulated light curves, both approaches of time delay 
determination give reasonable and compatible results. 



4. Results and Discussion 

The results for real data were obtained with the same procedure 
as was presented for the simulated light curves. Figure [5] shows 
the power spectrum P v calculated from the light curve of PKS 
1830-21 1. A periodic pattern is obvious on the distribution of P v 
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Fig. 3. Double power spectrum D a for the simulated lensed AGN Fig. 5. Measured power spectrum of PKS 1830-2 1 1 , plotted in 
of figure |2] D a is plotted in arbitrary units. The solid line is a fit arbitrary units, 
to a linear plus Gaussian profile. 
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Fig. 4. Autocorrelation function of the simulated lensed AGN of 
figure[2] The solid line is a fit to a linear plus Gaussian profile. 
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Fig. 6. Double power spectrum of PKS 1830-21 1 plotted in arbi- 
trary units. The solid line is a fit to a linear plus Gaussian profile. 



calculated from the PKS 1830-21 1 light curve. It is similar to the 
pattern expected from the simulations shown on figure [2] 

The autocorrelation function and the D a distribution calcu- 
lated for real data are shown on figures [7] and [6] A peak around 
27 days is seen in both distributions. Several other peaks are 
present on the autocorrel ation function as was already noted by 
Geiger & Schneider ( 1996) (see their Fig 1). The peak around 5 
days in the D a distribution is likely to be an artefact connected 
to the time variation of the exposure of the LAT instrument on 
PKS 1830-211. Using the method described in section [331 the 
significance of the peak around 27 days is found to be 1.1 cr in 
the autocorrelation function and 3 cr in the double power spec- 
trum D a . Fitting the position of the peak gives the time delay 
of a — 27.5 +1.3 days for the D a distribution. The fit of the 
autocorrelation function to a gaussian peak over an exponential 
background gives a = 27.07 ± 0.45 days. 

The double power spectrum distribution obtained for PKS 
1830-211 provides the first evidence for gravitational lensing 
phenomena in high energy gamma rays. The evidence is still at 
the 3 cr level but will likely improve by a factor of 2 over the 
lifetime of the Fermi satellite. Thanks to the uniform light curve 
sampling provided by Fermi LAT instrument, it is not necessary 
to identify features on the light curve to apply Fourier transform 
methods. The example of PKS 1830-21 1 shows that the method 
works in spite of the low photon statistic. Possible extensions of 
the present work are finding multiple delays in complicated lens 
systems or looking for unknown lensing systems in the Fermi 
catalog of AGNs. 




40 60 
Time lag [days] 



Fig. 7. Measured autocorrelation function of PKS 1830-21 1 . The 
solid line is a fit to an exponential plus Gaussian profile. 
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